#ifndef __MATCH_MIAN_CU__
#define __MATCH_MIAN_CU__

#include "GPU2BWT_kernel.h"
#include "SAResult.h"
#include "Match.cu"
#include "SAMatch.h"

#define MIN_NUM_SA (5)


typedef void (*KERNEL_FUNC)(uint3*, char* , int*, const int,  const unsigned int*, const int*,
                                        const unsigned int, BWT*, BWT *);


__device__ __host__
__inline__ bool smallSA(const int start, const int end) {
	return (end - start + 1) < MIN_NUM_SA ? (true) : (false);
}

__device__
__inline__ int addSA(uint3* result, char* saLabel,
			int* numResult, 
			const int maxNumOutput,
			const unsigned int queryId, 
			const unsigned int temp_start, 
			const unsigned int temp_end, 
			const char label) {
	if((*numResult) < maxNumOutput) {
		result[*numResult].x = queryId;
		result[*numResult].y = temp_start;
		result[*numResult].z = temp_end;
		saLabel[*numResult] = label;
		(*numResult) += 1;
	} else {//2. if the buffer is full
		*numResult = -1;
	}
	return *numResult;
}
	
__device__
bool new_bwtCalc(BWT *bwt, char baseId, 
		unsigned int * start, unsigned int * end) {
	if ((*start) > bwt->textLength || (*end) > bwt->textLength)
		return false;

        *start = GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue, (*start), baseId,bwt->inverseSa0) + 1;
        *end   = GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue,((*end)+1), baseId,bwt->inverseSa0);

					
	if (*start <= *end) {
		return true;
	} else {
		return false;
	}
}


__device__
bool new_ForwardSearch(BWT* revbwt, const unsigned int * query, 
			unsigned int lower, unsigned int upper, 
			unsigned int *temp_rev_start, unsigned int *temp_rev_end, 
			unsigned int *temp_start, unsigned int* temp_end,
			unsigned int *shareOccCount, 
			char* saOffset) {
	char c;
	unsigned int *occCount = shareOccCount;
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;

	*saOffset = -1;
	for (i=lower; i<upper && *temp_start<= *temp_end; i++) {

		c = getCharFromQueries(query, i);
		bwtCalcAll(revbwt, temp_rev_start, temp_rev_end, occCount_start, occCount_end, occCount);
		bwtCalcWithOccValue(revbwt->countValue, c, occCount_start[c], occCount_end[c], temp_rev_start, temp_rev_end);
		*temp_end = *temp_end - occCount[c];
		*temp_start = *temp_end - (*temp_rev_end - *temp_rev_start);

		if((*temp_start > 0) && (*temp_start <= *temp_end) && (smallSA(*temp_start, *temp_end))) {
			(*saOffset) = i;
			return true;
		}
	}
	if (*temp_start <= *temp_end && *temp_start > 0)
		return true;
	return false;
}


__device__ bool new_BackwardSearch(BWT *bwt, const unsigned int * query, 
				unsigned int lower, unsigned int upper, 
				unsigned int * temp_start, unsigned int * temp_end, 
				char* saOffset) {
	int i;
	char c;
	if(*temp_start > *temp_end)
		return false;

	for (i= (int)(upper-1); i>= (int)lower; i--) {
		c = getCharFromQueries(query, i);
		//if not found, then no results
		if(!new_bwtCalc(bwt, c,temp_start, temp_end)) {
			return false;
		}
	
		*saOffset = -1;
		//if the SA range is small enough, then no further 2bwt search
		//but you still should return true
		if(smallSA(*temp_start, *temp_end)) {
			*saOffset = i;
			return true;
		}
	}

	return true;
}


__device__ void new_BackwardSearch1Error(uint3* result, char* saLabel, 
		int &numResult, const int maxNumOutput,
		BWT * bwt, const unsigned int* theQuery,
		unsigned int queryId,
		unsigned int lower, unsigned int upper, 
		unsigned int * sa_start, unsigned int * sa_end) {
	unsigned int temp_start=(*sa_start),temp_end=*sa_end;
	int i;
	char c, misBase;

	//the buffer is exceeded
	if(numResult == -1) {
		return;
	}

	//otherwise
	for (i=upper-1; i >=(int)lower && *sa_start<=*sa_end; i--) {
		c = getCharFromQueries(theQuery, i);
		for(misBase=0; misBase < ALPHABET_SIZE; misBase++) {		
			if (misBase == c)
				continue;
			temp_start = *sa_start;
			temp_end = *sa_end;

			if(!new_bwtCalc(bwt, misBase, &temp_start, &temp_end)) {
				continue;
			} else {
				
				if(smallSA(temp_start, temp_end)) {
					//store the SA, but it is small SA
					if(-1 == addSA(result, saLabel, 
							&numResult, maxNumOutput, 
							queryId, temp_start, temp_end, 
							i)) {
						return;
					}
					continue; //no need for the further 2bwt for this misBase
				}
			}

			char saOffset = -1;
			if(new_BackwardSearch(bwt, theQuery, 0, i, 
					&temp_start, &temp_end, &saOffset)) {
			
				//store the SA, no matter it is a completed SA or not
				if(-1 == addSA(result, saLabel, &numResult, maxNumOutput, 
						queryId, temp_start, temp_end, 
						saOffset)) {
					return;
				}
			}
		}

		if(!new_bwtCalc(bwt, c, sa_start, sa_end)) {
			break;
		} else {
			if(smallSA(*sa_start, *sa_end)) {
				//store the SA
				if(-1 == addSA(result, saLabel, &numResult, maxNumOutput, 
						queryId, *sa_start, *sa_end, i)) {
					return;
				}
				break; // no need further SA
			}
		}
	}
}



__device__ int new_BWT1MismatchBackward(uint3* result, char* saLabel, 
		const int maxNumOutput,
		const unsigned int * query, const unsigned int queryId,
		BWT *bwt, BWT *revbwt) {
	unsigned int start= 1, end = bwt->textLength-1;
	char c;
	unsigned int i;
	unsigned long long packedPattern = 0;
	int queryLen = QUERY_LENGTH;

	for (i = 0; i < TABLE_SIZE; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, queryLen - TABLE_SIZE +i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];

	int numResult = 0;
	char saOffset;
	if(start <= end) {
		if(new_BackwardSearch(bwt, query, (QUERY_LENGTH/2), 
				queryLen-TABLE_SIZE, &start, &end, &saOffset)) {
		
			numResult = 0;
			if(smallSA(start, end)) {
				return addSA(result, saLabel, &numResult, maxNumOutput, 
						queryId, start, end, saOffset);
			}

			new_BackwardSearch1Error(result, saLabel, 
						numResult, maxNumOutput,
						bwt, query, queryId, 
						0, (QUERY_LENGTH/2),
						&start, &end);
		}	

	}

	return numResult;
}



__global__ void new_1MismatchBackward_kernel(uint3* d_result, char* d_saLabel,
		int* d_numOutputPerQuery, 
	 	const int maxNumOutputPerQuery,
		const unsigned int * Queries, const int* d_QueriesId,
		const unsigned int queryNum, BWT *bwt, BWT *revbwt) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;

	for(int idx = globalThreadOffset; idx < queryNum; idx += numTotalThread) {
		d_numOutputPerQuery[idx] = new_BWT1MismatchBackward(d_result + maxNumOutputPerQuery*idx,
							d_saLabel + maxNumOutputPerQuery*idx,
							maxNumOutputPerQuery, 
							Queries + QUERY_SIZE_IN_WORD*idx, 
							d_QueriesId[idx], 
							ptr_bwt, ptr_revbwt);
	}
}




__device__ void new_ForwardSearch1Error(uint3* result, char* saLabel, 
			int &numResult, const int maxNumOutput,
			BWT * revbwt, const unsigned int * query, unsigned int queryId,
			unsigned int lower, unsigned int upper,
			unsigned int * sa_rev_start, unsigned int * sa_rev_end, 
			unsigned int * sa_start, unsigned int * sa_end, 
			unsigned int * sharedOccCount) {
	unsigned int temp_start=*sa_start,temp_end=*sa_end;
	unsigned int temp_rev_start=*sa_rev_start, temp_rev_end=*sa_rev_end;
	unsigned int *occCount = sharedOccCount;
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;
	char c, misBase;
	
	if(numResult == -1) {
		return;
	}


	for (i = lower; (i < upper) && (*sa_start <= *sa_end); i++) {

		bwtCalcAll(revbwt,sa_rev_start, sa_rev_end,
				occCount_start, occCount_end, occCount);
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++) {
			if (misBase == c)
				continue;

			temp_rev_start = *sa_rev_start;
			temp_rev_end = *sa_rev_end;

			bwtCalcWithOccValue(revbwt->countValue, misBase,
					occCount_start[misBase], occCount_end[misBase], 
					&temp_rev_start, &temp_rev_end);

			temp_end = *sa_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start > temp_end)
				continue;

			//TODO: incorrect result?
			if(smallSA(temp_start, temp_end)) {
				if(-1 == addSA(result, saLabel, &numResult, maxNumOutput,
						queryId, temp_start, temp_end, i)) {
					return;
				}
				continue;
			}


			char saOffset = -1;
			if(new_ForwardSearch(revbwt, query, i+1,upper, &temp_rev_start, 
						&temp_rev_end, &temp_start,&temp_end,
						sharedOccCount + ALPHABET_SIZE, &saOffset)) {
				if(-1 == addSA(result, saLabel, &numResult, maxNumOutput,
						queryId, temp_start, temp_end, saOffset)) {
					return;
				}	
			}

		}

		bwtCalcWithOccValue(revbwt->countValue, c,occCount_start[c], occCount_end[c], sa_rev_start, sa_rev_end);
		*sa_end = *sa_end - occCount[c];
		*sa_start = *sa_end - (*sa_rev_end - *sa_rev_start);

		//TODO: incorrect results
		if((*sa_start <= *sa_end) && (smallSA(*sa_start, *sa_end))) {
			if(-1 == addSA(result, saLabel, &numResult, maxNumOutput,
					queryId, *sa_start, *sa_end, i)) {
				return;
			}
			break;
		}
	}
}


__device__ int new_BWT1MismatchForward(uint3* result, char* saLabel,
			const int maxNumOutput,
			const unsigned int * query, unsigned int queryId,
			BWT *bwt, BWT *revbwt, 
			unsigned int * sharedOccCount) {

	unsigned int start= 1, end = bwt->textLength-1, revStart, revEnd;
	char c;
	unsigned int i;

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= ( c & 3 );
	}
	r_packedPattern = l_packedPattern;
	l_packedPattern <<= (TABLE_SIZE - revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (TABLE_SIZE -revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(TABLE_SIZE - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	if (end - start < revEnd - revStart) {
		revStart =revEnd - (end - start );
	} else if (end - start > revEnd - revStart) {
		start = end - (revEnd - revStart);
	}

	char saOffset = -1;
	int numResult = 0;
	new_ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize, 
			(QUERY_LENGTH / 2), &revStart, &revEnd, 
			&start, &end,
			sharedOccCount, &saOffset);

	if(start <= end){
		if(smallSA(start, end)) {
			return addSA(result, saLabel, &numResult, maxNumOutput, 
						queryId, start, end, saOffset);
		}

		new_ForwardSearch1Error(result, saLabel, 
					numResult, maxNumOutput,
					revbwt, query, queryId, 
					(QUERY_LENGTH/2), QUERY_LENGTH,
					&revStart, &revEnd, &start, &end, 
					sharedOccCount);
	}

	return numResult;
}



__global__ void new_1MismatchForward_kernel(uint3* d_result, char* d_saLabel, 
			int* d_numOutputPerQuery, 
			const int maxNumOutputPerQuery,
			const unsigned int * Queries, const int* d_QueriesId,
			const unsigned int queryNum, 
			BWT *bwt, BWT *revbwt) {
	extern __shared__ unsigned int sharedOccCount[];
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;

	for(int idx = globalThreadOffset; idx < queryNum; idx += numTotalThread) {
		d_numOutputPerQuery[idx] = new_BWT1MismatchForward(d_result + maxNumOutputPerQuery*idx,
								d_saLabel + maxNumOutputPerQuery*idx,
								maxNumOutputPerQuery,
								Queries + QUERY_SIZE_IN_WORD*idx, 
								d_QueriesId[idx], 
								ptr_bwt, ptr_revbwt,
								sharedOccCount + OCC_COUNT_PER_THREAD * threadIdx.x);
	}
}



__device__ int new_2MismatchCaseA(uint3* result, char* saLabel, 
			int maxNumOutput,
			const unsigned int * query, unsigned int queryId, 
			BWT *bwt, BWT *revbwt) {
	unsigned int temp_start=0,temp_end=0;
	unsigned long long packedPattern = 0;
	unsigned int i;
	char c, misBase;
	unsigned int start, end;
	
	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, QUERY_LENGTH - TABLE_SIZE + i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];


	int numResult = 0;

	if(start <= end &&  BackwardSearch(bwt, query, sizeY, QUERY_LENGTH - TABLE_SIZE, &start, &end)) {
		for (i=sizeY-1; i>=1 && start<=end; i--) {
			c = getCharFromQueries(query, i);
			for(misBase = 0; misBase < ALPHABET_SIZE; misBase++) {
				if( misBase == c)
					continue;
				temp_start = start;
				temp_end = end;
				if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){
					new_BackwardSearch1Error(result, saLabel,
							numResult, maxNumOutput - numResult,
							bwt, query, queryId,
							0,i,
							&temp_start, &temp_end);
					//the buffer is full
					if(numResult == -1) {
						return -1;
					}				
				}
			}
			bwtCalc(bwt, c, &start, &end);
		}

	}

	return numResult;
}



__global__ void new_2MismatchCaseA_kernel(uint3* d_result, char* d_saLabel,
				int* d_numOutputPerQuery, const int maxNumOutputPerQuery,
				const unsigned int * d_Queries, const int* d_QueriesId,
				const unsigned int queryNum, 
				BWT *d_bwt, BWT *d_revbwt) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = globalThreadOffset; idx < queryNum; idx += numTotalThread) {

		d_numOutputPerQuery[idx] = new_2MismatchCaseA(d_result + idx*maxNumOutputPerQuery,
								d_saLabel + idx*maxNumOutputPerQuery,
								maxNumOutputPerQuery,
								d_Queries + QUERY_SIZE_IN_WORD*idx, 
								d_QueriesId[idx], 
								&bwt, &revbwt);
	}
}


__device__ int new_2MismatchCaseB(uint3* result, char* saLabel, 
				int maxNumOutput,
				const unsigned int * query, unsigned int queryId, 
				BWT *bwt, BWT *revbwt, 
				unsigned int * sharedOccCount) {
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
	unsigned int *occCount = sharedOccCount;
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	//unsigned long long mask;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;
	int numOutput = 0;

	l_packedPattern = 0;
	r_packedPattern = 0;
	rev_packedPattern = 0;
	// look-up the first characters
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= (c & 3);
	}
	r_packedPattern = l_packedPattern;
	// if the look-up tables are of different size
	l_packedPattern <<= (bwt->lookupTable->tableSize -  revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end  = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end - start<revEnd-revStart) {
		revStart = revEnd - (end - start);
	} else if (end - start > revEnd - revStart) {
		start = end - (revEnd - revStart);
	}

	if (start <= end && ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize,sizeY,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE)) {
		for (i=sizeY; i<sizeZ && start<=end; i++) {
			bwtCalcAll(revbwt, &revStart, &revEnd, occCount_start, occCount_end, occCount);
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++) {
				if(misBase == c)
					continue;
				temp_rev_start = revStart;
				temp_rev_end = revEnd;

				bwtCalcWithOccValue(revbwt->countValue,misBase,
						occCount_start[misBase], occCount_end[misBase], 
						&temp_rev_start, &temp_rev_end);

				temp_end = end - occCount[misBase];
				temp_start = temp_end - (temp_rev_end - temp_rev_start);

				if(temp_start <= temp_end){

					new_ForwardSearch1Error(result, saLabel, 
							numOutput, maxNumOutput - numOutput,
							revbwt, query, queryId, 
							i + 1, sizeZ, 
							&temp_rev_start, &temp_rev_end,
							&temp_start, &temp_end,
							sharedOccCount + ALPHABET_SIZE);
					if(numOutput == -1) {
						return -1;
					}
				}
			}

			bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

			end = end - occCount[c];
			start = end - (revEnd - revStart);
		}
	}

	return numOutput;
}


__global__ void new_2MismatchCaseB_kernel(uint3* d_result, char* d_saLabel, 
					int* d_numOutputPerQuery, const int maxNumOutputPerQuery,
					const unsigned int * d_Queries, const int* d_QueriesId,
					const unsigned int queryNum, 
					BWT *d_bwt, BWT *d_revbwt) {
	extern __shared__ unsigned int sharedOccCount[];
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = globalThreadOffset; idx < queryNum; idx += numTotalThread) {
		d_numOutputPerQuery[idx] = new_2MismatchCaseB(d_result + idx*maxNumOutputPerQuery,
								d_saLabel + idx*maxNumOutputPerQuery, 
								maxNumOutputPerQuery,
								d_Queries + QUERY_SIZE_IN_WORD*idx, 
								d_QueriesId[idx], 
								&bwt, &revbwt,
								sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);
	}
}


/**
collect the SA result and unfinished queries
TODO: most jobs should be done on the GPU directly
@return the number of undone reads
*/
int collectResultAndQuery(SAResult &saResult,
			SAMatch &saMatch,
			const int numWordPerQuery,
			unsigned int* d_outQuery, 
			int* d_outQueryId,
			const unsigned int* d_inQuery,
			const int* d_inQueryId,
			const int numInQuery,
			const uint3* d_result, 
			int* d_numSAPerQuery, 
			const char* d_saLabel,
			const int maxNumSAPerQuery) {

	//get the number of SA
	int numSA = 0;
	int* h_numSAPerQuery = (int*)malloc(sizeof(int)*numInQuery);
	FROMGPU(h_numSAPerQuery, d_numSAPerQuery, sizeof(int)*numInQuery);
	for(int i = 0; i < numInQuery; i++) {
		if(h_numSAPerQuery[i] != -1) {
			numSA += h_numSAPerQuery[i];
		}
	}

	//collect the correct SA and small SA
	uint3* h_outSA = saResult.getSAForAppend();
	uint3* h_sa = (uint3*)malloc(sizeof(uint3)*maxNumSAPerQuery*numInQuery);
	FROMGPU(h_sa, d_result, sizeof(uint3)*maxNumSAPerQuery*numInQuery);
	char* h_saLabel = (char*)malloc(sizeof(char)*maxNumSAPerQuery*numInQuery);
	FROMGPU(h_saLabel, d_saLabel, sizeof(char)*maxNumSAPerQuery*numInQuery);
	int count = 0;
	int countSmallSA = 0;
	for(int i = 0; i < numInQuery; i++) {
		for(int j = 0; j < h_numSAPerQuery[i]; j++) {
			const int k = i*maxNumSAPerQuery + j;
			if(h_saLabel[k] == -1) {
				h_outSA[count] = h_sa[k];		//the correct, complete SA
				count += 1;
			} else if(h_saLabel[k] >= 0) {
				saMatch.add(h_sa[k], h_saLabel[k]);	//the small SA
				countSmallSA += 1;
			} else {
				ERROR_EXIT("unknown saLabel");
			}
		}
	}
	saResult.numSA += count;
	free(h_sa);
	free(h_saLabel);
	
	cout << "correct SA: " << count << endl;
	cout << "Small SA  : " << countSmallSA << endl;

	//collect the unfinished queries
	int numOutQuery = 0;
	unsigned int* h_inQuery = new unsigned int[numWordPerQuery*numInQuery];
	FROMGPU(h_inQuery, d_inQuery, sizeof(unsigned int)*numWordPerQuery*numInQuery);
	int* h_inQueryId = new int[numInQuery];
	FROMGPU(h_inQueryId, d_inQueryId, sizeof(int)*numInQuery);
	for(int i = 0; i < numInQuery; i++) {
		if(h_numSAPerQuery[i] == -1) {
			numOutQuery += 1;
		}
	}
	unsigned int* h_outQuery = new unsigned int[numWordPerQuery*numOutQuery];
	int* h_outQueryId = new int[numOutQuery];
	count = 0;
	for(int i = 0; i < numInQuery; i++) {
		if(h_numSAPerQuery[i] == -1) {
			h_outQueryId[count] = h_inQueryId[i];
			memcpy(h_outQuery + count*numWordPerQuery, h_inQuery + i*numWordPerQuery, 
					sizeof(unsigned int)*numWordPerQuery);
			count++;
		}
	}
	TOGPU(d_outQuery, h_outQuery, sizeof(unsigned int)*numWordPerQuery*numOutQuery);
	TOGPU(d_outQueryId, h_outQueryId, sizeof(int)*numOutQuery);


	free(h_numSAPerQuery);
	delete[] h_inQuery;
	delete[] h_inQueryId;
	delete[] h_outQuery;
	delete[] h_outQueryId;

	return numOutQuery;
}


void cuda_mismatch(SAResult &saResult, SAMatch &saMatch,
			QueryHelper& helper, BWT* d_bwt, BWT* d_revbwt, 
			MATCH_CASE matchCase) {

	KERNEL_FUNC kernelFunc;
	if(matchCase == MIS1_BACKWARD) {
		cout << "===> MIS1_BACKWARD" << endl;
		kernelFunc = new_1MismatchBackward_kernel;
	} else if(matchCase == MIS1_FORWARD) {
		cout << "===> MIS1_FORWARD" << endl;
		kernelFunc = new_1MismatchForward_kernel;
	} else if(matchCase == MIS2_CASE_A) {
		kernelFunc = new_2MismatchCaseA_kernel;
	} else if(matchCase == MIS2_CASE_B) {
		kernelFunc = new_2MismatchCaseB_kernel;
	} else {
		ERROR_EXIT("NOT SUPPORTED YET!");
	}

	int maxNumSAPerQuery = 2; //(helper.GPUResultBufSize)/(helper.numQuery)/sizeof(uint3);
	int numQuery = helper.numQuery;

	//prepare the memory
	cutilSafeCall(cudaMemcpy(helper.d_queryBuf, helper.d_query, 
				sizeof(unsigned int)*helper.numWordPerQuery*helper.numQuery, 
				cudaMemcpyDeviceToDevice));
	cutilSafeCall(cudaMemcpy(helper.d_queryIdBuf, helper.d_queryId, 
				sizeof(int)*helper.numQuery, 
				cudaMemcpyDeviceToDevice));
	unsigned int* d_inQuery = helper.d_queryBuf;
	int* d_inQueryId = helper.d_queryIdBuf;	
	unsigned int* d_outQuery = helper.d_queryBuf_temp;
	int* d_outQueryId = helper.d_queryIdBuf_temp;
	const int numBlock = 512;
	const int numThread = 128;
	const int sharedMemSize = numThread * OCC_COUNT_PER_THREAD * sizeof(int);

	unsigned int timer = 0;

	int passId = 0;
	while(numQuery > 0) {
		cout << "pass " << passId << " .............." << endl;
		cout << "numQuery = " << numQuery << endl;
		cout << "maxNumSAPerQuery = " << maxNumSAPerQuery << endl;

		//kernel
		timer = 0;	
		startTimer(&timer);
		kernelFunc<<<numBlock, numThread, sharedMemSize>>>(
					helper.d_result, helper.d_saLabel,
					helper.d_numSAPerQuery, 
					maxNumSAPerQuery,
					d_inQuery, d_inQueryId,
					numQuery, d_bwt, d_revbwt);
		cutilCheckMsg("kernelFunc");
		cutilSafeCall(cudaThreadSynchronize());
		endTimer(&timer, "2BWT search kernel");

		//collect the SA results and unfinished queryies
		timer = 0;
		startTimer(&timer);
		numQuery = collectResultAndQuery(saResult,
				saMatch,
				helper.numWordPerQuery, 
				d_outQuery, d_outQueryId,
				d_inQuery, d_inQueryId,
				numQuery,
				helper.d_result, helper.d_numSAPerQuery, 
				helper.d_saLabel,
				maxNumSAPerQuery);
		cudaThreadSynchronize();
		endTimer(&timer, "collectResultAndQuery");

		//update the control parameters for the next pass
		swap(d_inQuery, d_outQuery);	
		swap(d_inQueryId, d_outQueryId);
		maxNumSAPerQuery *= 4;

		passId += 1;

		cout << endl;
	}
}


#endif /*__MATCH_MIAN_CU__*/
